Resolving all-order method convergence problems for atomic physics applications 
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The development of the relativistic aU-order method where all single, double, and partial triple 
excitations of the Dirac-Hartree-Fock wave function are included to all orders of perturbation theory 
led to many important results for study of fundamental symmetries, development of atomic clocks, 
ultracold atom physics, and others, as well as provided recommended values of many atomic prop- 
erties critically evaluated for their accuracy for large number of monovalent systems. This approach 
requires iterative solutions of the linearized coupled-cluster equations leading to convergence issues 
in some cases where correlation corrections are particularly large or lead to an oscillating pattern. 
Moreover, these issues also lead to similar problems in the CI-|-all-order method for many-particle 
systems. In this work, we have resolved most of the known convergence problems by applying 
two different convergence stabilizer methods, reduced linear equation (RLE) and direct inversion of 
iterative subspace (DIIS). Examples are presented for B, Al, Zn^, and Yb"''. Solving these conver- 
gence problems greatly expands the number of atomic species that can be treated with the all-order 
methods and is anticipated to facilitate many interesting future applications. 

PACS numbers: 31.15.bw, 31.15.ac, 06.30.Ft, 31.15.ag 
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I. INTRODUCTION 



The coupled-cluster (CC) method has been success- 
fully applied to solve quantum many-body problems in 
quantum chemistry P, |3] as well as computational 
atomic 0j and nuclear physics A relativistic lin- 

earized variant of the coupled-cluster method (which is 
numerically symmetric and is generally referred to as 
"all-order method") was developed for atomic physics 
applications in Refs. 043 ■ I* '^^^ of the most accu- 
rate methods currently being used in the atomic struc- 
ture calculations. This approach was extremely success- 
ful and led to accurate predictions of energies, transi- 
tion amplitudes, hyperfine constants, polarizabilities, C3 
and Cq coefficients, isotope shifts, and other properties 
of monovalent atoms, as well as the calculation of parity- 
nonconserving (PNC) amplitudes in Cs, Fr, and Ra+ (see 
[sl-fisj and references therein). Further development of 
the all-order approach, that included triple excitations 
and non-linear terms yielded the most precise evaluation 
of the PNC amplitude in Cs [isj and consequent re- 
analysis of Cs experiment [16| . This work provided the 
most accurate low-energy test of the electroweak sector of 
the Standard Model to date, placed constraints on a va- 
riety of new physics scenarios beyond the SM, and, when 
combined with the results of high-energy collider experi- 
ments, confirmed the energy dependence (or "running") 
of the electroweak force over an energy range spanning 
four orders of magnitude (from '-^lO MeV to ~100 GeV). 
AU-order method was also used for development of ultra- 
precise atomic clocks (l7l-2lL ultracold atom and quan- 
tum information studies 2 2142 q and many other appli- 
cations. We refer the reader to review [13] for details 
of the all-order method and its applications. The all- 



order method is also used as a part of the Cl-f all-order 
approach for study of more complicated systems [27| . 



The all-order method requires iterative solutions of the 
linearized coupled-cluster equations leading to conver- 
gence issues in some cases when correlation corrections 
are very large or produce an oscillating iterative pattern. 
The initial guess of the solution is based on the low-order 
perturbation theory. Therefore, if high-order correlation 
corrections are large, initial guess is very poor leading 
to very slow convergence or failure of the straightforward 
iterative scheme. In addition, initial non-linear CC equa- 
tions may have more than one solution, so a convergence 
to non-physical solutions may occur. Several such prob- 
lems have been identified over the years and led to failure 
to apply all-order approach for many important applica- 
tions. For example, all or almost all of the low-lying nd 
and nf states of B, Al, Zn+, Cd+, Hg+, and Yb+ do not 
converge if standard Jacobi-type iterative procedure is 
applied. In the case of Yb+, even core equations do not 
converge. Convergence problems also cause complete fail- 
ure of the all-order approach for super-heavy elements, 
such as element 113 (eka-Tl). All these convergence is- 
sues in monovalent systems lead to the same problems in 
the application of the Cl+all-order approach [27] to the 
corresponding divalent systems, such as Al"*", Hg, Yb, 
etc. since this method required prior solution of LCCSD 
equations for one-particle orbitals. There are several in- 
teresting present applications of these atoms and ions 
that require high-precision calculations possible with all- 
order techniques. For example, several of these systems 
are used or proposed for optical clocks [13, [28l - [3l| requir- 
ing precise knowledge of the blackbody radiation (BBR) 
shift which is hard to accurately measure. BBR shift is 
a leading source of uncertainties for many of the atomic 
clock schemes. Yb is used for an ongoing PNC experi- 
men t 13211 as well as studies of degenerate quantum gases 
[33I |34| owing to a number of available isotopes. The 
best available Yb PNC amplitude value is only accurate 
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to 20% m. 

The convergence issues that arise in the solutions of 
eigenvalue equations have a long history in general quan- 
tum chemistry and several methods have been developed 
to address them [36.- 42J . Most of these methods are based 
on the fundamental idea of effective reduction of original 
large functional space and solution of the projected to the 
reduced (Krylov) subspace of the simplified equations. 
This idea was implemented for the first time in a quan- 
tum chemical application by Lanczos 4^] , who facilitated 
a partial diagonalization of a large matrix by transform- 
ing to a much smaller Krylov subspace, followed by a 
matrix triangularization procedure. In the present work, 
we consider two such convergence techniques, namely, re- 
duced linear equation (RLE) [sB, |33] and the direct in- 
version of iterative subspace (DIIS) [H, [3^. Both of 
the methods use approximate solutions obtained from 
few last iterations as Krylov reduced functional subspace 
onto which the linearized equations are projected and in 
which the projected system of equations is solved. The 
convergence of the methods is based on the construction 
of error vectors. Different choices of the error vectors 
lead to different implementation of the methods. Among 
the most popular error vectors are: 1) the difference of 
subsequent iterations and 2) "true" error vector (e.g. dif- 
ference between exact solution and it's approximation). 
In our work, both the convergence methods use the same 
best least squares approximation to the true error vector 
and thus are rather relative. Moreover, our variant of 
DIIS can be regarded as a "symmetric" version of RLE 
(see below). However, while DIIS method is chosen to 
minimize the error vector in the least-squares sense, the 
RLE differs from it by requiring that this vector within 
the basis vanishes. We formulate here implementations 
of the RLE and DIIS methods for our variant of the 
coupled-cluster equations and test these stabilizer meth- 
ods on several specific examples, in which we were able 
to resolve the convergence problems listed above. We 
also studied the effectiveness of these two techniques in 
solving specific types of the convergence problems as well 
as accelerating convergence in all other cases. Accelera- 
tion of convergence is particulary important for further 
CI-|-all-order use since it requires solving all-order equa- 
tions for a large number of one-particle orbitals. 

Below, we briefiy outline the essence of the convergence 
stabilization procedures. In the coupled-cluster method, 
the desired exact wave function jV') is obtained by ap- 
plying (a yet unknown operator) exp(T) on some refer- 
ence wave function |(/)), for example, the Dirac-Hartree- 
Fock (DHF) wavefunction. For a closed-shell system 
with N electrons, the cluster operator T = ^Tp (where 
p = 1, 2, 3..., N) has the form: 

Tp = l/p\ ^ Pmn...a6...aLaJj...aaaf,... (1) 
mn..ah.. 

Here, orbitals m,n... are single-particle excited states; 
a, 6... are core states which are occupied in p's are 
cluster amplitudes (also called excitation coefficients) 



FIG. 1: (Color online) The failure of the LCCSD straight- 
forward iteration procedure for the 3s state in boron. The 
calculated correlation energy is plotted as a function of the 
iteration number. The dashed (red) line indicates the value 
of the experimental correlation energy. 



and a* and a are creation and annihilation operators 
with respect to the quasi- vacuum state Finally, p 
is the number of core electrons excited when applying 
Tp to In the linearized coupled-cluster single-double 
(LCCSD) method, only Ti and T2 are retained and non- 
linear terms in the expansion of exp(T) are truncated. 
The LCCSD equations are conventionally solved by an it- 
erative scheme, symbolically written as p("+i) = F(p("^), 
F being specified later in Section |TT1 In this paper, this 
type of straightforward iteration procedure is referred to 
as the conventional iterations scheme (CIS). 

Both RLE and DIIS convergence stabilization proce- 
dures form solution as the linear combination of 
cluster amplitudes (p("\ p^"'"-^^ p^""')) accumulated 
from I previous CIS iterations. Further details of the 
LCCSD method and RLE and DIIS schemes are dis- 
cussed in Sections HIl and [Till 

An example of failed conventional iteration procedure 
is shown in Fig. [U where we plot the LCCSD correla- 
tion energy, SE, as a function of a number of valence 
LCCSD iterations for the 3s state of boron. The experi- 
mental correlation energy (—0.0079754 a.u.) is indicated 
by the horizontal dashed line. It is computed by sub- 
tracting DHF energy from the experimental result. The 
LCCSD 3s correlation energy diverges from the experi- 
mental values dramatically and begins to oscillate after 
a number of iterations. The convergence criteria is set 
to terminate the iteration procedure when the relative 
difference between two consecutive iterations is reduced 
below 0.00001. The convergence is not reached even after 
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500 iterations. As demonstrated below, this problem is 
completely resolved by the use of either RLE or DIIS pro- 
cedures and convergence to the above criteria is reached 
within 30 iterations. 

This paper is organized as follows: in Section we 
describe the LCCSD method and the conventional iter- 
ation procedure (CIS) of solving the LCCSD equations. 
In Section Hill we formulate RLE and DIIS schemes for 
LCCSD equations. In Section IIVI we analyze perfor- 
mance of the RLE and DIIS schemes for various cases. 
Finally, in Section V, we draw the conclusions. 



II. LINEARIZED SD COUPLED-CLUSTER 
METHOD 

In the present implementation of the CC method, the 
exact valence wave function is obtained from the 
lowest-order DHF state, 

\<P.)=al%), (2) 

by applying a wave operator J7 = Af[exp(r)] 

|^,) = ri|0„), (3) 

where |0c} is the core DHF state and N[..] designates the 
normal product of operators with respect to a closed-shell 
core. Taking into account only the Ti and T2 terms in 
Eq. (dJ, and truncating i7 past the linear terms in the 
expansion of the exponential leads to the LCCSD ansatz 
for the wave operator 



2 ^ 

■'Tia mnab 



(4) 



Here, Sc and Dc {Sy, Dy) are the core (valence) single 
and double terms, respectively. 

To find the cluster amplitudes (or excitation coefB- 
cients) p, we need to specify the Hamiltonian. In our 
approach, we use the Hamiltonian \2} H — Hq + G: 



H = 



'^eiN[alai] + gijkiN[ala]aiak] , (5) 



ijkl 



where Hq is the one-electron lowest-order DHF Hamilto- 
nian and G is the residual Coulomb interaction. Indices i, 
j, k, and I range over all possible single-particle orbitals, 
and gijki are the two-body Coulomb matrix elements. A 
set of coupled equations for the cluster operators (r)„: 

(Tc)i = 5c, (r„)i = Sy, (Tc)2 = Dc, and (Ty)2 = 

may be found from the Bloch equation Q. For monova- 
lent systems [43 [: 

(e„ — Ho){Tc)n 
(e„ -I- SEy — Ho){Ty)n 



{QGil} 

connected, n i (6) 

{QGil} 

connected, n 7 (7) 



where SE^ = {(f)y\Gil\(j)v) is the valence correlation en- 
ergy and Q — 1 — is the projection operator. 
Note that Eq. ([6|) contains only the core cluster opera- 
tors, while Eq. ([7]) contains both core and valence cluster 
operators. The core equations ([6]) are solved first, and 
the resulting CC core amplitudes are subsequently frozen 
and used in the valence equations ([7]). 

The summations over the magnetic quantum numbers 
TO in Eqs. ((H) and ([7]) are performed analytically. After 
the angular reduction, the equation for the reduced single 
core cluster amplitudes p{ma) takes form [§, [i^ : 

(ca - eyn)p{ma) = (8) 

no V 

(_l)ia+Jb+ic+J„ 

- > > r—TTTi Zk{cbna)pk{nmcb) 

k neb ^■^''J^^J 

+ 2^2^ r— |7u {mbrn)pk {rnab) } . 

k rnh L nj L J 

Here, [j] = 2j + 1, k is the relativistic angular momen- 
tum quantum number, p(ma) and pk{mnab) are reduced 
single and double cluster amplitudes, Xk{mnab) are re- 
duced two-body Coulomb matrix elements, and 



Zk{mnab) 



Xk{mnab) 

jm ja ^ 



The equations for the reduced double core cluster am- 
plitudes pk{mnab) are given by: 



{£ab - £m«) Pkimnab) = Xk{mnab) 
+ AiXi(cdab)pk' {mncd) 

cd l.k' 

+ ^ ^ A2X1 {■mnrs)pk' (rsab) 



(9) 



rs Lk' 





a-i^b 


+ 


m n 



^Xk{mnrb)p(ra)S^^^^ +^Xk{cnab)p{mc)S^^^^ 



(_l)jc+J,.+fe 

y , j^jj Zk{cnrb)pk{mrac) 



where Ai are angular coefhcients given in j45| and = 
Ei+Ej. The valence equations have exactly the same form 
as the core equations with the replacement of index a by 
the valence index v everywhere and an addition of the 
valence correlation energy SEy into the energy difference 
on the left-hand side, i.e., {Sa — Sm)pi'ma) — > (£,„ — e,„ 4- 
SEy)p(rnv). 

Implementation of the RLE and DIIS procedures re- 
quires rewriting the equations for the cluster amplitudes 
in a specific vector form. We introduce the vector nota- 
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tion: 



p(ma) 
Pk{mnab) 



where p{ma) and pk{mnab) are to be understood as 
columns composed of all amplitudes for single and double 
excitations, respectively, i.e., for all possible values of m, 
n, o, 6, and k indexes. Then, the core equations given by 
Eqs. ([5]) and ^ may be combined as 



D t= a A t 



(10) 



where 



a = 





Xk (mnab) 



D 



and A • t includes all terms on the right-hand sides of the 
Eqs. ^ and (O except for Xk{mnah), which is included 
in a. 

Valence equations may be written in the same way with 

p{mv) 
Pkimnvb) 



and 





Xk{mnvb) 



D 



The main difference between the core and valence equa- 
tions for the implementation of the RLE and DIIS is the 
dependence of the valence array D on the iteration num- 
ber, since 5Ev is recalculated after every iteration. In the 
core case, D remains constant. 
Solving Eq. pO)) for t gives 



-D"\a + A-t) 



(11) 



The above equation can be solved iteratively as 

t(™+i) = -D-\a+ A-t(")). (12) 

The iteration usually starts by inserting t'"-* =0 on the 
right hand side of Eq. and finding t'-^-'. As we demon- 
strated in Fig.[Tl convergence of this straightforward iter- 
ative scheme is occasionally very slow or fails altogether. 
The convergence methods that we develop in the next 
section will alleviate such problems and lead to faster 
convergence rates. 

III. RLE AND DIIS 



next best solution of Eq. p2|) . The new answer is then 
used for another initialization of the CIS and the two 
steps are repeated until convergence to specified criteria 
is reached. In this section, we present the general RLE 
and DIIS formulas and derive their explicit form for the 
LCCSD equations. 

After accumulating m + 1 iteratively found solu- 
tions, t^^-', t^^-*, t^™"''^-', next best approximation can 
be found as their linear combination, 



I -1-1] 



a it 



(13) 



The quantities ai are the weights that have to be de- 
termined by solving a system of equations constructed 
from previously found m-\-\ CIS solutions. We note that 
t(™+i) is not included in the linear combination (|13p . 
but is used to find ai coefficients. Therefore, we use the 
notation t'™"'"-'^! instead of t'-™"'"^' to distinguish between 
the m -f- 1*'* solution found through the use of RLE/DIIS 
methods and the initial CIS result, respectively. 

Both direct inversion of iterative space (DIIS) and 
reduced linear equation (RLE) methods seek to mini- 
mize the error between the iteratively found solutions 
of Eq. (Ilip and the exact answer. The error minimiza- 
tion is the basis for finding the appropriate Ui to form 
the approximate solution tf'"+^l. Both methods also use 
a least square approach to the error minimization. Since 
the exact answer is unknown, approximations are used in 
the minimization process. The approximate solution, as 
mentioned before, is constructed as a linear combination 
of a series of iteratively found solutions. The difference 
between the DIIS and the RLE methods is in the assump- 
tions they make in order to minimize the errors. Further 
details of the difference between the two methods and 
derivations of the DIIS/RLE formulas can be found in 
the Appendix A. 

We rewrite Eq. (dU]) as 

a+ (A + D)t = a + Bt = 0. (14) 
The DIIS formula for determining cr, is given by Eq. (jA7l) : 

T'^B'^a + T^B^BTcr = . (15) 
The RLE formula for determining (7^ is given by Eq. (|A9I) : 
T^(a + BTcr) = 0. (16) 



Both Eqs. p5)) and (|T6|) can be written as a system of m 
equations: 



In this section, we formulate implementation of RLE 
and DIIS methods for the LCCSD equations (ITT]) dis- 
cussed in the previous section. Both methods are two- 
step procedures. In the first step, a few iterative solutions 
t*^'^ of Eq. (IT^ are found (same as the CIS). In the second 
step, a linear combination of these t^*) is used to find the 



a + Rcr = . (17) 

Solving the above system of equations for a can be easily 
done with standard linear algebra methods. The result- 
ing coefficients (Ji are substituted into Eq. to obtain 
best new approximate solution tt™+^l. 
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Next, we write R and a of Eq. p7|) in their explicit 
forms for DIIS and RLE methods. Substituting A — D 
for B into DIIS equation ([T5|) yields for the cji 

t''«(A + -Dfa + t''«(A + D)'^(A + D)t«a,; = . 



Using Eq. we find that A • t(*) = -(D • t^'+i) + 

a). Replacing the dot products involving A with ones 
involving D yields explicit form of DIIS matrix for core 
orbitals 



k ' "fe "k 



k 

"k '-fe 



(18) 



The RLE equations for core orbitals are obtained by 
repeating the same steps as for the DIIS approach but 
starting from Eq.([T6]). The resulting RLE equations for 
R and a are 



% = Y.*k i^ki + Dki)t 



kl 



k 



(i+i) 

k 



akt'i^] 



(19) 



We noted in the previous section that D depends on 
the correlation energy, 5Ey^ in the case of the valence 
equations leading to the dependence of D on the iteration 
number. Therefore, the substitution D — D^*^ must 
be made to rewrite the DIIS and RLE equations above 
for the valence orbitals. To derive the final form of the 
equations, we have to introduce a somewhat arbitrary dot 
product and normalization definitions. The explicit form 
of the core RLE equations is obtained by substituting the 
expressions for D, a, and t from the previous section into 
Eq. (dH): 

ma 

+ E X! ^^''b ~ P'^L (mnab) 

L mnab 

X [p^' (mnab) — p^'l^^^ (mnab)] — ai , 
ai = — E] E/ TTT"^-^ (mnab) p^l^ [mnab) . (20) 
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FIG. 2: (Color online) Comparison of the RLE5, DIIS5, and 
DIIS7 schemes for the 3s state of boron. The correlation 
energy is given in a.u. 

RLE equations for valence case are given by 



ai = 



ma 

X [p(j') (ma) - (ma)] 

+ X! XI Mfl if"'' ~ ^ ^^v^) P^L {mnab) 

L mnab J L 

X [p^i^ (mnab) — p^i^^^ {mnab)] — ai , 
E E 1^^^ (mna6) p^'' {mnab) , 



(21) 



L mnab 



where [L] =2L + 1. 

The implementation of the RLE and DIIS methods 
proceeds as follows. In step one, our code makes a limited 
number, m + 1, of LCCSD iterations using the CIS. This 
is done to find m+1 series of single and double cluster am- 
plitudes, p^^\ma) and p^j^\mnab) that are then saved. In 
step two, a separate subroutine applies the DIIS or RLE 
equations to these stored cluster amplitudes to find the 
appropriate R and a matrices. The m-dimensional lin- 
ear equation pT|) is solved for a. The next best solution 
of the LCCSD equations is then found by substituting a 
into Eq. These two steps are repeated until con- 

vergence is reached according to a specified criteria. In 
the next section, we discuss the results of the application 
of the DIIS and RLE procedures to the solution of the 
LCCSD equations in the cases that do not converge or 
converge to non-physical answers with the conventional 
iteration scheme. 



IV. RESULTS AND DISCUSSION 

In this section, we study and compare the convergence 
characteristics of the DIIS and the RLE methods. We 
include a number of test cases in four different systems. 
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TABLE I: Convergence tests of the LCCSD equations with 
CIS, RLE, and DIIS methods for B and Al. CIS is the con- 
ventional iterations scheme (no convergence stabilizer). RLE5 
designates RLE convergence method with 5 pre-stored itera- 
tions. Last column gives resulting correlation energy in a.u. 
*Cases where maximum number of iterations allowed during 
run was reached. 



Atom 


State 


Method 


# of iter. 


Converged? 


5Ev{a.u.) 


B 


Core 


CIS 


21 


Yes 








RLE5 


9 


Yes 








DIIS5 


15 


Yes 




B 


2pi/2 


CIS 


18 


Yes 


-0.0293907 




RLE5 


13 


Yes 


-0.0293906 






DIIS5 


13 


Yes 


-0.0293908 


B 


3s 


CIS 


70* 


No 


-0.0091643 






RLE4 


70* 


No 


-0.0070438 






DIIS4 


70* 


No 


-0.0089454 






RLE5 


30 


Yes 


-0.0089491 






DIIS5 


22 


Yes 


-0.0089472 






DIIS7 


31 


Yes 


-0.0089488 


B 


3pi/2 


CIS 


44 


Yes 


-0.0056292 




RLE5 


23 


Yes 


-0.0056294 








19 


Yes 


A r\r\ ca^o A 






RLE8 




Yes 








Dims 




1 eo 


-u.uutJUZryo 


B 


3^3/2 


CIS 


85* 


No 


f\ f\ClC^ A A cir\ 

-0.0884489 




mm? 


71 


Yes 


u.uou t ooo 






RLE8 


UU 


Yes 


U.UUU ( tJOU 






-L/-L100 


o ± 


Yes 




A 1 

Al 


Core 


CIS 


13 


Yes 








RLE5 


8 
O 


Yes 




Al 


3pi/2 


CIS 


16 


Yes 


-0.0245810 




RLE5 


11 


Yes 


-0.0245798 






DIIS5 


14 


Yes 


-0.0245811 


Al 


4s 


CIS 


20 


Yes 


-0.0079907 






RLE5 


13 


Yes 


-0.0079906 






DIIS5 


18 


Yes 


-0.0079905 


Al 


3^3/2 


CIS 


70* 


No 


-0.0209573 




DIIS6 


89 


Yes 


-0.0208637 






RLE8 


86 


Yes 


-0.0208662 






DIIS8 


49 


Yes 


-0.0208662 


Al 


4^3/2 


CIS 


70* 


No 


-0.0213727 




RLE8 


300* 


No 


0.0022280 






DIIS8 


81 


Yes 


0.0022478 






DIIS9 


91 


Yes 


0.0022477 


Al 


4d5/2 


CIS 


70* 


No 


-1.3775005 




RLE8 


165 


Yes 


0.0022737 






DIIS8 


97 


Yes 


0.0022710 






DIIS9 


73 


Yes 


0.0022712 


B, Al, 


Zn+, 


and Yb+ 


which 


aave a large number of 



states that do not converge with the conventional iter- 
ative scheme (CIS). We also test the ability of the RLE 
and the DIIS to accelerate convergence in the cases where 
the CIS does converge. The main purpose of these tests 
is to provide general guidelines of how to accelerate or to 
achieve convergence using the RLE and the DIIS meth- 



ods. The conclusions and observations of this section 
may be extrapolated to other systems for both all-order 
and CI-|~all-order approaches. 

The summary of B and Al convergence tests is given in 
Table H] We find that convergence patterns for two fine- 
structure multiplet states, for example 3pi/2 and 3p3/2 
states, are generally very similar. Therefore, we list only 
TT'Pi/2 and states with the exception of the Ad states 

of Al. Tests were performed for both states of the mul- 
tiplet as an additional check, since similar results are 
expected. The results are given for the '2,pi/2, 3s, ?>Pi/2, 
and 3c?3/2 states of B and the 3pi/2, 4s, 3^3/2, 4c?3/2 
and 4(^5/2 states of Al. The resulting LCCSD correla- 
tion energy is listed in the last column of the table in 
a.u. The convergence method is specified in the third col- 
umn. CIS refers to the initial straightforward iteration 
scheme. RLE5 designates the RLE convergence method 
with 5 pre-stored CIS iterations. Similarly, DIIS8 refers 
to the DIIS convergence scheme with 8 pre-stored itera- 
tions. The fourth column indicates the iteration number 
at the end of the run. Cases where the maximum num- 
ber of iterations allowed during the run was reached are 
marked with asterisk. In these cases, convergence did not 
occur. The convergence criteria was set to terminate the 
iteration procedure when the relative difference between 
two consecutive iterations is reduced below 0.00001. The 
same convergence criteria is used in all valence test runs. 
Only the core and the 2p states of B converge with the 
CIS. In the case of Al, all nd states do not converge with 
the CIS. All of the cases in Table U converge with the 
DIIS8. We note that we did not hst the RLE5 and the 
DIIS5 results for many of the nd states because conver- 
gence was not achieved. In cases where all methods lead 
to convergence, both the RLE and the DIIS have accel- 
erated convergence rates relative to the CIS. 



0.001 




5 10 15 20 25 30 35 
Number of Iterations 

FIG. 3: (Color online) Comparison of the DIIS5, RLE8, and 
DIIS8 schemes for the 3ci3/2 state of boron. The correlation 
energy is given in a.u. 

We may draw two general conclusions from our tests: 
1. If a particular LCCSD run converges with the CIS, 
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then the RLE5 appears to be the most efficient 
method in accelerating the convergence. 

2. If a particular LCCSD run does not converge with 
the CIS, the DIIS8 or the DIIS9 appear to be the 
most efficient in attaining and accelerating the con- 
vergence. 

We note that these two rules are not absolute, but they 
serve to be good initial guidelines. Our further tests on 
other (much heavier) systems confirmed these guidelines. 
We note that the RLE5 is not sufficient to achieve conver- 
gence for most of the divergent cases. The only exception 
in Table U is the 3s state of B. However, while the CIS 
never converges to our standard criteria for the 3s state, 
it nearly converges to correct result before exhibiting di- 
verging and oscillating pattern of Fig. [TJ In this case, 
accumulation of only 5 iterations is sufficient. However, 
in the case of the nd states, the CIS is never close to con- 
verging to a correct number and subsequently the RLE5 
does not work. Occasionally, DIIS9 may achieve conver- 
gence where DIIS8 would not. Using even larger number 
of stored iterations does not improve convergence or ef- 
ficiency. DIIS10-DIIS12 runs for the Ad states converged 
to non-physical answers in two instances, but to correct 
results in all other cases. Number of iterations varied 
significantly from case to case. The results of all con- 
verged runs listed in Table |T] are consistent within the 
convergence criteria, as expected. 

We implemented the DIIS/RLE strategies for two sep- 
arately developed LCCSD codes. The calculations were 
carried out usin g tw o different finite basis sets, the B- 
splines of Ref. [4a| and the dual-kinetic-basis sets of 
Ref. [13] ■ Even though the basis sets and the conver- 
gence criteria used for each code made slight differences 
in the values, the general observations on the convergence 
patterns remains the same. 

We illustrate different convergence patterns of the RLE 
and the DIIS methods for the 3s and 3^3/2 states of boron 
in Figs. [5] and [31 In Fig. [21 the values of the correla- 
tion energies for the 3s states of boron are obtained from 
different schemes that are listed on the graph. RLE5, 
DIIS5, and DIIS7 results after = 20 interactions are 
indistinguishable at this plot scale and are not shown. 
These schemes converge after 30, 22, and 31 iterations, 
respectively. While the CIS results appear close to con- 
verged value, convergence was never reached and corre- 
lation energy began to oscillate as illustrated in Fig. [H 
The RLE5 and the DIIS5 results are identical to the CIS 
ones for the first four iterations. The 5th value is differ- 
ent for three of the schemes as this (m -I- l)-th value (see 
Eq. (fT3)) ') is replaced by the RLE or the DIIS predictions 
for RLE5 and DIIS5. We observe that these predictions 
are significantly closer to converged result than the 5th 
CIS iteration. After that, the RLE5 and DIIS5 results 
are sharply adjusted at TV = 10 when the second call to 
the RLE/DIIS stabilizer codes is made. The DIIS7 be- 
havior is similar to the one just described, except that 
it accumulates 7 iterations before the DIIS procedure is 



TABLE II: Convergence tests of the LCCSD equations with 
CIS, RLE, and DIIS methods for Zn+ and Yb+. CIS is 
the conventional iterative scheme (no convergence stabilizer) . 
DIIS8 designates the DIIS convergence method with 8 pre- 
stored iterations. Last column gives resulting correlation en- 
ergy in a.u. *Cases where maximum number of iterations 
allowed during run was reached. 



Atom 


State 


Method 


# of iter. 


Converged? 


5Ey{a.u.) 


Zn+ 


5pi/2 


CIS 


39 


Yes 


-0.0066119 




uiiby 




Yes 


-U.UUDDUoo 


Zn+ 


4^3/2 


CIS 


70* 


No 


-0.0977215 




RLE5 


67 


Yes 


-0.0045508 






PlTTQQ 

Ulloo 




Yes 


-U.UU4ooii 


Zn+ 


4^5/2 


CIS 


70* 


No 


-0.1149058 




RLE5 


93 


Yes 


-0.0045266 








1 


Yes 


-0. 004026 1 


Zn+ 


5^3/2 


CIS 


70* 


No 


-0.1936330 




RLE5 


28 


Yes 


-0.0018929 






DIIS8 


18 


Yes 


-0.0018942 


Zn+ 


4/5/2 


CIS 


200* 


No 


-0.0008697 




DIIS8 


200* 


No 


-0.0007949 






DIIS9 


153 


Yes 


-0.0007948 


Zn+ 


4/7/2 


CIS 


70* 


No 


-0.0007970 




DIIS8 


173 


Yes 


-0.0007891 






DIIS9 


195 


Yes 


-0.0007891 


Yb+ 


Core 


CIS 




No 








RLE5 


12 


Yes 








DIIS5 


12 


Yes 





invoked and now the 7th value gets much closer to the 
final answer. 

In Fig.[3l the values of the correlation energies obtained 
from CIS, DIIS5, RLE8, and DIIS8 are plotted for the 
3(i3/2 states of boron. The RLE8 and the DIIS8 results 
after = 35 appear identical on the graph at this scale 
and are not shown. The RLE8 and the DIIS8 converge to 
our criteria after 66 and 51 iterations, respectively. Very 
similar behavior of the RLE8 and the DIIS8 is observed, 
with the RLE8 energy oscillations being slightly larger 
after the RLE subroutine pass. However, other tests show 
that the RLE8 in general converges slower, sometimes 
dramatically so, than the DIIS8. The CIS values diverge 
completely and increase rapidly. The DIIS5 seems to be 
converging at TV = 35, but does not in fact reach selected 
criteria even after 100 iterations. 

The summary of the selected Zn"*" and Yb^ conver- 
gence tests is presented in Table [Hi The results are given 
for the 5pi/2, 4^3/2, 4^5/2, and 5^3/2, 4/5/2, and 4/7/2 
states of Zn+ and the Yb+ core. The 4s and Apj states of 
Zn+ and the 6s, 6pj, 7s, and 5dj states of Yb"*" converge 
with the CIS, so we have omitted these results from the 
table. However, it is worth pointing out that RLE5 ac- 
celerates convergence for all these states compared to the 
CIS. Table [Til demonstrates that the DIIS reduces num- 
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TABLE III: Comparison of B, Al, Zn+, and Yb^ removal energies (in cm ^) with experiment [j^ ]. Rows labeled "Dif." give 
relative difference with experimental values in %. 



B 


2pi/2 


2p3/2 


3s 


3pi/2 


3p3/2 


3^3/2 


34/2 




Expt. 


-66928 


-66913 


-26888 


-18316 


-18314 


-12160 


-12160 




SD 


-67049 


-67035 


-27105 


-18497 


-18495 


-12494 


-12494 




Dif. 


-0.18% 


-0.18% 


-0.81% 


-0.99% 


-0.99% 


-2.7% 


-2.7% 




Al 




4s 


"J ""3/2 


'-'"'5/2 


Tj 

^^3/2 


^""5/2 






rjxpL . 


-^oz / o 


-zzyoi 














SD 


-48271 


-23069 


-17295 


-17289 


-6652 


-6647 






Dif. 


0.02% 


-0.60% 


-8.4% 


-8.4% 


-9.1% 


-9.1% 






Zn+ 


4s 






'-'f 1/2 


Ari', In 
^^Aj 1 


^""5/2 


^ 


^7 7/2 












( you 


47009 


9fiQ1 ^ 


-Z ( uuu 




1 AARRA 

- L'±^\JO'± 


QR991 






47Q9Q 
-"4: ( yzy 


-4: 1 OOU 




-Z ( uoo 












47QQ4 


47Q4fi 


o«Q9Q 


-z ( uoo 


Dif. (SD) 


0.14% 


0.20% 


0.19% 


0.24% 


0.11% 


0.11% 


0.09% 


-0.02% 


Dif. (SDpT) 


-0.23% 


-0.15% 


-0.15% 


0.08% 


-0.02% 


-0.03% 


-0.02% 


-0.01% 


Yb+ 


6s 


6pi/2 


6p3/2 


7s 


5^3/2 


5(^5/2 


5/5/2 


5/7/2 


Expt. 


-98207 


-71145 


-67815 


-43903 


-75246 


-73874 


-27704 


-27627 


SD 


-98961 


-71016 


-67480 


-44060 


-76141 


-74700 


-28080 


-28062 


SDpT 


-99107 


-71084 


-67592 


-44115 


-77764 


-76317 






Dif. (SD) 


-0.77% 


0.18% 


0.49% 


-0.36% 


-1.19% 


-1.12% 


-1.36% 


-1.57% 


Dif. (SDpT) 


-0.91% 


0.09% 


0.33% 


-0.48% 


-3.2% 


-3.2% 








■7 I I I I I I I I I I I 

' 1 2 3 4 5 6 7 8 9 10 11 12 

Number of Iterations 



FIG. 4: (Color online) Comparison of the CIS, RLE5, and 
DIIS5 schemes for the Yb^ core. The correlation energy is 
given in a.u. RLE5 and DIIS5 data appear identical at this 
scale and are shown as a single curve. 



ber of iterations for the 5pi/2 states by a factor of 3 or 
better. Zn+ and Yb"*" tests confirm our conclusions (1) 
and (2), on the previous page. We were unable to achieve 
convergence for higher l-pj states in Yb"*". This problem 
is not present in Zn+, where LCCSD for the 5p states 
converges even with CIS as shown in Table |TT1 Perhaps 
other convergence approaches are needed to resolve this 
issue. 

The case of Yb^ core is particularly interesting, since 



core iterations generally converge well with the CIS. Yb+ 
core is an exception, however, most likely due to very 
large 4/ shell contributions that lead to oscillation of 
the correlation energy. We plot the CIS, the RLE5, and 
the DIIS5 results for the Yb+ core correlation energy in 
Fig. m The RLE5 and the DIIS5 results appear identical 
at this scale and are shown as a single curve. Both meth- 
ods are successful at fixing the CIS's oscillation problem. 

The comparison of the B, Al, Zn+, and Yb+ removal 
energies with experiment [48[ is given in Table Hill Rows 
labeled "Dif." give relative difference with experimental 
values in %. The energies here are given in cm~^ . Most of 
these states did not converge with the CIS, so it is impor- 
tant to establish the accuracy of this approach for such 
cases. Breit interactions and contributions from higher 
partial waves are also included. The B and Al ioniza- 
tion potentials, B 2p3/2, and Al 4s SD energies are in 
agreement with experiment. We consider only monova- 
lent states for all of these systems. The SD approxima- 
tion does not account for mixing with the core-excited 
states such as 3s3p^ in Al. Therefore, larger disagree- 
ment with experiment is expected in cases where mixing 
with these hole-two-particle states is large. A particular 
example is the 3d and 4c? states of Al. The lower 3s^nd 
levels heavily mix with the 3s3p^ levels. However, 
the mixing coefficient for this configuration never exceeds 
30%. As a result, these levels are distributed over sev- 
eral lower nd levels, resulting in two sets of levels being 
listed as Ss^M 48, 49] {[y^D] and pD]). In TableHIIl 
we compare the 4d results with the second sets of levels 
i?D]). 

We also included partial valence triples perturbatively 
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(LCCSDpT) to investigate if the LCCSDpT method 
would improve the theory-experiment agreement for Zn+ 
and Yb+ . This method is described in detail in [9] . Since 
triple equations are not explicitly iterated in this ap- 
proach, implementation of the RLE and the DIIS method 
is exactly the same as in the SD code. Convergence tests 
of the LCCSDpT method exhibit essentially the same 
pattern as the tests of the LCCSD method discussed 
above, and a similar number of iterations was generally 
required for LCCSD and LCCSDpT calculations for the 
same states run with the same parameters. 

As shown in Table IIIH we find an excellent agreement 
of all Zn"*" data with experiment. Inclusion of perturba- 
tive triples somewhat improves the agreement with ex- 
periment for most states. The accuracy decreases for 
Yb^, as expected, owing to much softer and heavier core 
and strong mixing of monovalent states with one-hole- 
two-particle states in this system. Nevertheless, for Yb"'" 
the average accuracy for removal energies is at the level 
of 1% (see Table mi) 



V. CONCLUSION 

We have successfully implemented the RLE and DIIS 
convergence techniques in the LCCSD and LCCSDpT 
methods for high-precision atomic many-body calcula- 
tions. Most of the convergence problems were resolved 
using these methods. Acceleration of convergence was 
demonstrated for all cases where all-order equations con- 
verge with straightforward iteration scheme. Numerous 
tests were performed to establish general recommenda- 
tions for the RLE/DIIS use for various purposes. We 
find that if particular case converges with CIS, RLE5 
appears to be the most efficient in achieving and accel- 
erating convergence. If particular case does not converge 
with CIS, DIIS8 or DIIS9 appear to be the most efficient 
in accelerating convergence. Solving these convergence 
problems greatly expands the number of atomic species 
that can be treated with the all-order methods and is 
anticipated to facilitate many interesting future applica- 
tions for studies of fundamental symmetries as well as 
atomic clock and ultracold atom research. 



Appendix A 

Derivations of general formulas in this Appendix 
mainly follows Appendix of Ref. [36] . Consider solving a 
general linear equation of the form 

a + Bt = 0, (Al) 

which is a system of linear equations of dimension k with 
vector t being the exact solution that we would like to 
find. We make the best approximation to the exact so- 
lution by using m(< k) nonorthogonal and linearly inde- 
pendent vectors T = (t^^^t^^), ...,t(")), where each t^') 



is a fc-dimensional vector. We find this best approxima- 
tion as a linear combination of t'''''s: 

m 

^ J2 C7,t(') = 0- • T . (A2) 

Here, ai are the weights of the optimized solution that 
needs to be determined. We note that t'''"^^' is not in- 
cluded in the linear combination (IA2p . but is used to 
construct matrices such as shown in Eas. (|^D)) and (|2ip . 
Therefore, t'™+^'' needs to be also found through the CIS. 
Therefore, we use the notation t['"+-^l instead of 1^™+^) to 
distinguish between the m + solution found through 
the use of RLE/DIIS methods and the CIS result, respec- 
tively. 

First, we try to derive an ideal equation to find UiS 
as if we know the exact solution to Eq. (jAip . To find 
the best approximation, we need to minimize the error 
between the approximate and the exact answers. To this 
end, we use the least square optimization approach. The 
error is e = t — t["'+-'^l. The least square optimization of 
E = e^e with respect to a then yields 

^ = -2T^(t-T-a) =0. (A3) 
After solving for a and substituting it in Eq. (|A2p , we get 
^[m+i] ^ T(T^T)-iT'^t . (A4) 

However, not knowing what the exact solution t is, 
the above formula is of little use. The DIIS and RLE 
are based on replacing t with approximations. Substi- 
tuting tl^+^l instead of t in Eq.dJT]), will make Ea. (|Il|) 
inhomogeneous : 

a-HBt[™+il =a + BT-(T = e. (A5) 

where e is a vector with constant elements. 

The difference between the RLE and DIIS methods is 
in their choice of error to minimize, e. The DIIS takes 
the error to be e of Eq. I|A5|) . Then to get the best ap- 
proximation, we need to minimize E = e^e with respect 
to a: 

— ^2(-BT)^(a + BTa) ==0. (A6) 

OCT 

Therefore, the coefficients a that lead to the best approx- 
imation satisfy the DIIS equation; 

T'^B'^a + T^B^BTcr = . (A7) 

The RLE requires that the best least square approxi- 
mation e[™+il to e vanishes in the space of T. Following 
the structure of Eq. (jA4|) : 

g[m+i] ^ T(T^T)-iT^e = T(T^T)-iT^(a + BTa) = . 

(A8) 
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Since T is made of linearly independent vectors, el'^+^l 
is only zero if: 



T^(a + BT(T) =0. 



(A9) 



Eqs. and for DIIS and RLE, respectively, 

correspond to Eqs. (fT5|l and p6|) in the paper. 
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